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Abstract — The  paper  deals  with  state  estimation  of  nonlin¬ 
ear  non-Gaussian  systems  with  a  special  focus  on  the  Gaussian 
sum  filters.  To  achieve  a  higher  estimate  quality,  state  and 
measurement  predictive  moments  appearing  in  the  filters  are 
computed  by  the  randomized  unscented  transform,  which 
provides  asymptotically  exact  estimates  of  the  moments.  The 
use  of  the  Gaussian  sum  filter  employing  the  randomized 
unscented  transform  is  introduced  and  the  proposed  algo¬ 
rithm  is  illustrated  in  a  numerical  example.  The  analysis 
of  the  numerical  example  involves  a  comparison  of  several 
filters  using  a  number  of  performance  metrics  both  absolute 
and  relative,  assessing  the  point  estimate  quality,  the  estimate 
error  quality,  and  the  density  estimate  quality. 

Keywords:  state  estimation,  nonlinear  filtering,  non- 
Gaussian. 

I.  Introduction 

State  estimation  of  nonlinear  non-Gaussian  discrete¬ 
time  stochastic  dynamic  systems  is  of  utter  importance 
in  many  areas  such  as  target  tracking  [1],  [2],  navigation, 
signal  processing,  fault  detection,  and  adaptive  and  optimal 
control  problems. 

The  Bayesian  framework  is  a  suitable  tool  for  the  state 
estimation  problem  as  it  enables  a  non-Gaussianity  analysis 
with  respect  to  the  description  of  random  quantities.  A  gen¬ 
eral  solution  to  recursive  state  estimation  problems  within 
the  Bayesian  framework,  is  given  by  the  Bayesian  recursive 
relations  (BRR’s).  They  produce  probability  density  func¬ 
tions  (pdf’s)  of  the  state  conditioned  by  the  measurements. 
The  pdf’s  represent  a  full  stochastic  description  of  the  state, 
which  itself  is  unmeasurable. 

The  closed-form  solution  to  the  BRR’s  is  available 
only  for  a  few  special  cases,  e.g.,  for  a  linear  Gaussian 
system  [3],  where  the  BRR’s  solution  corresponds  to  the 
Kalman  filter.  When  the  measurement  model  or  the  state 
transition  model  is  nonlinear,  an  approximate  method  must 
be  applied.  These  approximate  methods  can  be  divided  into 
two  groups  with  respect  to  the  validity  of  the  resulting 
estimates  [4].  The  first  group  of  methods  provides  results 
with  validity  within  a  neighborhood  of  the  point  estimate 
only  and  thus  they  are  called  local  methods.  The  second 


group  of  methods  provides  results  valid  within  almost  the 
whole  state  space  and  are  called  global  methods. 

The  local  methods  provide  approximate  conditional 
mean  and  associated  covariance  matrix  of  the  estimate 
error.  Popular  approximations  to  nonlinear  functions  used 
by  the  methods  are  (i)  the  Taylor  series  expansion  of  the 
nonlinear  functions  in  the  system  description  leading  to  the 
extended  Kalman  filter  (EKF)  or  the  second-order  extended 
Kalman  filter  [5],  (ii)  the  polynomial  linearization  of  the 
nonlinear  function  by  a  first  or  second  order  polynomial 
interpolation  [6],  [7]  leading  to  the  divided  difference  filters 
(DDF’s),  (iii)  the  stochastic  linearization  [8]  approximating 
a  random  variable  by  a  set  of  points,  which  are  transformed 
through  nonlinear  functions  and  leading  to  the  unscented 
Kalman  filter  (UKF)  [9] — [1 1],  and  (iv)  calculating  approx¬ 
imate  moments  of  a  random  variable  by  means  of  the 
Gauss-Hermite  quadrature  [6]  or  cubature  integration  rules 
and  leading  to  quadrature  or  cubature  filters  [12],  In  [13] 
a  randomized  unscented  Kalman  filter  has  been  proposed 
(RUKF)  based  on  a  degree  3  stochastic  integration  rule 
(SIR3),  which  is  combination  of  the  cubature  rule  and  the 
Monte  Carlo  (MC)  method. 

Due  to  the  approximation  of  the  state  estimate  condi¬ 
tional  pdf  by  the  first  two  moments  only,  the  local  methods 
are  not  very  practical  for  non-Gaussian  problems.  Here, 
the  global  estimation  methods,  providing  an  approximation 
of  the  full  conditional  pdf,  achieve  a  good  quality  perfor¬ 
mance.  However,  note  that  the  improved  performance  is 
usually  paid  by  higher  computational  costs.  There  are  three 
main  approaches  to  the  global  filtering  method  design:  (i) 
the  analytical  approach  based  on  Gaussian  sum  approxima¬ 
tion  of  pdf’s  [14],  [15]  and  using  approximation  techniques 
of  the  local  methods,  (ii)  the  numerical  approach  using 
point-mass  approximation  of  the  conditional  state  pdf’s  and 
solving  the  integrals  in  the  BRR’s  numerically  [16],  [17], 
and  (iii)  the  simulation  approach  taking  advantage  of  the 
BRR’s  solution  by  the  MC  methods  and  approximating  the 
conditional  state  pdf  by  an  empirical  representation  [18], 

[19]. 

The  analytical  global  methods  are  based  on  a  multiple 
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application  of  local  methods,  e.g.  the  Gaussian  Sum  Filter 
(GSF)  consists  of  a  bank  of  EKF’s,  or  the  Sigma  Point 
Gaussian  Sum  Filter  [20]  using  a  bank  of  DDF’s  or  UKF’s. 
In  problems  involving  multiple  mode  pdf’s,  the  analytical 
global  methods  provide  a  high  estimate  quality  as  they  are 
able  to  follow  the  modes  with  individual  local  methods. 

The  goal  of  the  paper  is  to  propose  a  filter  based  on 
Gaussian  sum  representation  of  pdf’s  and  utilizing  the 
randomized  unscented  transform  (RUT)  which  has  been 
used  in  [13]  as  a  basis  for  designing  the  RUKF.  The  RUT 
will  be  used  to  compute  statistics  of  a  transformed  random 
variable.  Comparing  to  the  unscented  transform  (UT), 
appearing  in  the  UKF,  the  RUT  offers  asymptotically  exact 
estimates  of  the  statistics.  The  RUT  gives  a  flexible  and 
computationally  efficient  method  for  computing  posterior 
moments  when  measurement  and/or  kinematic  models  are 
nonlinear. 

The  paper  is  organized  as  follows:  System  specification 
and  nonlinear  non-Gaussian  state  estimation  by  means  of  a 
generic  Gaussian-sum  based  global  filter  are  introduced  in 
Section  II.  The  RUT  is  discussed  in  Section  III.  Section  IV 
deals  with  the  algorithm  of  the  proposed  filter  with  a 
special  focus  on  implementation  details.  In  Section  V,  the 
proposed  filter  will  be  applied  to  a  numerical  example  for 
performance  evaluation  and  comparison,  and  concluding 
remarks  are  drawn  in  Section  VI. 


II.  System  Specification  and  generic 
Gaussian-sum  based  global  filter 

Let  the  discrete-time  nonlinear  stochastic  system  be 
considered  in  the  following  form 

X£+1  =  f k(*k)  +  Wfc,  k  —  0, 1, 2, . . . ,  (1) 

z*  =  h*(x*)  +  v*,&  =  0, 1,2, (2) 


where  the  vectors  xk  e  M."x  and  zk  e  represent  the 
unmeasurable  state  of  the  system  and  measurement  at  time 
instant  k,  respectively,  fk  :  -*  R"x  and  — »■ 

M"1  are  known  vector  functions,  and  wk  e  R"x,  v*  e 
are  the  state  and  measurement  white  noises,  which  are 
supposed  to  be  mutually  independent.  The  pdf’s  of  the 
noises  p(wk)  and  p(vk)  are  assumed  to  be  known.  The 
initial  state  xo  is  independent  of  the  noises  and  its  pdf  is 
supposed  to  be  known. 

The  general  solution  to  the  estimation  problem  (i.e. 

finding  x*  based  on  knowledge  of  zk  =  [zo,  zi, . . . ,  z^]) 
is  given  by  the  BRR’s  computing  pdf’s  of  the  state  condi¬ 
tioned  by  the  measurements  [3].  These  pdf’s  provide  a  full 
description  of  the  estimated  state.  The  BRR’s  are  assumed 
in  the  following  form 


p(xk+l\zk) 


p(xk\zk  l)p( zk\xk) 
p(zk\zk-') 


j  p(xk+l\xk)p(xk\zk)dxk. 


(3) 

(4) 


where  ^(x^lz^)  is  the  filtering  pdf,  p(xk+  \\zk)  is  the  one- 
step  ahead  predictive  pdf,  p(xk- fi|x^)  and  p(zk\xk)  are  the 


state  transition  pdf  obtained  from  (1)  and  the  measurement 
pdf  obtained  from  (2),  respectively,  and  p(zk\zk~l)  — 
J  p(xk\zk~l)p(zk\xk)dxk.  The  closed  form  solution  to  the 
BRR’s  is  available  only  for  a  few  special  cases  [3],  In 
other  cases  it  is  necessary  to  apply  an  approximation  in 
the  BRR’s  solution. 

The  Gaussian  sum-based  global  filter  requires  specifi¬ 
cation  of  the  noises  and  initial  state  pdf’s  in  the  form  of 
Gaussian  sums  as  follows: 

N0 

p(xi o)  =  yiao}^{xo;  Xq^Po0}’  (5) 

i=l 

qk 

ptyk)  =  w['},  q[°],  (6) 

i=  1 
n 

=  (7) 

i= 1 

where  yV{y;  y,  Py]  denotes  Gaussian  distribution  of  a  ran¬ 
dom  variable  y  parametrized  by  its  mean  y  and  covariance 
matrix  Py.  The  parameters  ai'\  and  yk  are  positive 
weights  of  particular  Gaussian  terms  with  their  sum  being 
equal  to  1.  All  the  parameters  (i.e.  the  weights,  means  and 
covariance  matrices)  are  assumed  to  be  known. 

If  any  of  the  noises  or  the  initial  state  has  a  distribution 
different  from  Gaussian  sum  distribution,  its  Gaussian  sum 
approximation  has  to  be  found  (e.g.  using  the  expectation- 
maximization  (EM)  method  [21]).  It  has  been  proved  that 
the  approximation  by  a  Gaussian  sum  can  be  arbitrarily 
accurate  [3],  Considering  the  noises  given  by  the  Gaussian 
sum  pdf  (5-7),  the  following  algorithm  of  a  global  filter 
provides  a  generic  solution  to  the  estimation  problem. 

Algorithm  1:  Generic  Gaussian  Sum  Global  Filter 

Step  1:  ( initialization )  Set  the  time  instant  k  =  0  and 
define  a  priori  initial  condition  p(xo|z_1)  =  p(x o)  as  a 
sum  of  (Vo|-i  =  No  Gaussian  terms. 

Step  2:  (filtering )  The  filtering  pdf  is  approximated  by 

Nk\k 

p(xk\zk)  &  '^a^kJ\f{xk;  x[lMkl),  (8) 

i=  1 


where  Nk\k  —  Nk\k-\  ■  rk.  The  filtering  mean  x[‘[.  and  the 

covariance  matrix  P^J.  of  the  ;-th  term  J\f{xk;x^k,P^k} 
are  computed  using 


_  £(/)  4-K(''(7;  —7^  1 

xk\ k  ~  xk\k-l  +  ^klk^k  £k\k-l>’ 

(9) 

p(0  pO)  lfCOpO)  (tt (0  \T 

rk\k  ~  rk\k-\  ^k^zMk-t^  k\k>  ’ 

(10) 

where  =  P(^k]k_ ,  (p %k_ , )  and 

^k\k-\  — 

(11) 

pzV.  =  covS-,[h^)]+Rf> 

(12) 

Ptlk\k-1  =  Ek\k-iUxk  ~  *lt|ifc-i)(hfc(xfc)  - 

-t)T]> 

(13) 
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where  E^j<_]  [■]  is  a  shorthand  notation  for 


Again,  the  indices  j  and  /  are  given  by 


E^(x,;x«  ^ l}[']  =  /  [^{x^-pP^-lidXA. 


™d  cov”’_,[  ]  for  COV^.,0,^ 

The  likelihood  of  the  /th  term  with  respect  to  the 
known  measurement  zk  is  given  by 

p!Vi}-  a4) 

The  indices  j  and  /  are  given  by 


j  —  i  —  L— - \Nm-X,  j  =  1, . . . ,  Nk\k-i,  (15) 

Nk\k-i 

/  =  1  +  L77— ^-J,  l  =  l,...,rk  (16) 

Nk\k- 1 

and  i  =  l,...,Nk\k-  The  symbol  [x]  denotes  the  floor 
function,  i.e.  the  largest  integer  less  than  or  equal  to  x. 
The  filtering  weight  «[',!.  of  each  term  is  given  by 


j  —  i  ~  L— — \Nk\k,  j  —  1,  •  •  ■ ,  Nk\k, 

Nk\k 

i  —  1 

i  —  1  +  L— — J,  l  =  i,---,qk 
Nk\k 

and  i  —  1, . . . ,  Nk+i\k- 

Let  k  —  k  +  1,  Then  go  to  Step  2. 


III.  Randomized  Unscented  Transform 

The  RUKF  proposed  in  [13]  employs  the  RUT  as  a  mean 
to  calculate  approximate  predictive  statistics  (state  means 
and  covariance  matrices)  of  the  state  and  measurement.  The 
RUT  is  a  special  case  of  SIR3  proposed  in  [24],  [25],  The 
SIR3  aims  at  evaluating  an  integral  of  the  form 

p  =  f  tp(x)  (2jt)-^/2  |  P|“  1/2e“  2  (x~*)7p~'  A-^dx, 


a(‘)  =  _ ak\k-Jk  Ck _ 

Uk\k  Nt[k-1  VI1  (j)  (l)ANk\k-l(J-l)+0  '  K  > 

2Lj= 1  Zw=l  ak\k-Uk  <=k 

Step  3:  ( global  point  estimate )  The  global  filtering  mean 
and  covariance  matrix  can  be  obtained  by  the  following 
relations 

Nk\k 

*k\k  =  E[x/:|zA']  =  (18) 

(=i 

p  k\k  =  covfelz*] 

JVitl* 

—  [p[|1  +  (xfc|i-x[p)(x*p-x[|].)rj  .  (19) 

i=l 

Similarly,  other  point  estimates  (e.g.  a  mode  or  median) 
can  be  obtained  from  the  global  state  estimate  (8). 

Step  4:  ( reduction )  Generally,  the  state  or  measurement 
noise  with  the  Gaussian  sum  distribution  causes  an  expo¬ 
nential  growth  of  the  number  of  Gaussian  terms  in  the  sum 
(8),  which  must  be  reduced  to  keep  computational  costs 
reasonable  [22],  [23], 

Step  5:  ( prediction )  The  predictive  pdf  is  approximated  by 

Nk+l\k 

p(xk+l\zk)  «  ^  0EJfc+t|Jfc'^’[X*:-t-i;  Xjt-f-ilifc*  Pfc+il*}’  (20) 
;= 1 

where  Nk+i\k  =  Nk\k  •  qk  and  o^+ip  =  1  The 

particular  predictive  mean  x^^  and  covariance  matrix 

Pfe+tl*  are  comPuted  by 

*1+11*  =  1  Cxfc)]  +  w[\  (21) 

pS1|*  =  covS-i[f*(*t)]  +  Q*)-  (22) 


where  <p(-)  is  an  arbitrary  function.  Note  that  relation  (25) 
can  be  interpreted  as  computation  of  the  expected  value 
of  the  function  <p(x)  where  x  is  a  random  variable  with 
p(x)  —  «Af{x;  x,  P],  i.e. 

=  £',Ar{x;x,P}[<P(x)].  (26) 

The  SIR3  proposed  in  [24]  for  a  solution  to  (25)  is  given 
by  the  following  algorithm: 


sjorithm 


stochastic  integration  rule 


Step  1:  Choose  an  error  tolerance  e  and  a  maximum 
number  of  iterations  Nmax. 

Step  2:  Initialize  the  number  of  iterations  N  —  0,  initial 
value  of  the  integral  p  =  0„r  x  1 ,  and  initial  square-error 
of  the  integral  estimate  V  =  0„vX„r  and  calculate  a  square 
root  VF7  of  P(  such  that  Pv  =  VPvVPv"1-  Note  that  0„x/, 
denotes  a  a  x  b  matrix  of  zeros. 


Step  3:  Repeat  (until  N  —  Nmax  or  ||V||  <  e)  the 

following  loop: 


a)  Set  N  =  N  +  1. 

b)  Generate  a  uniformly  random  orthogonal  matrix 
U  of  dimension  nx  x  nx  and  generate  a  random 
number  p  from  chi-squared  distribution  p  ~ 
Chi(«A  +  2). 

c)  Compute  a  set  of  points  X;  and  appropriate 
weights  a>i  according  to 


X.o  —  x,  coo  —  1  *2 , 

X,  =  X  -  P\Ky/Vx)u  COi  =  ^2, 

lknx+i  —  X  +  p\i(\fPx)i,  COnx-\-i  —  COj, 


where  i  —  1,2 ,...,nx  and  the  term  (yl*v)- 
represents  the  / -th  column  of  the  matrix  VFv- 
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Figure  1.  Comparison  of  points  and  weights  generated  by  the  CR 
and  SIR3  for  M2  state  space:  left  -  comparison  of  the  weights,  right  - 
comparison  of  the  points  placement  (dots  -  SIR3,  stars  -  CR). 

d)  Compute  the  value  S  of  the  integral  at  current 
iteration 

2nx 

S  =  2><P(X,), 

1=0 

and  use  it  to  update  the  approximate  value  p  of  the 
integral  and  square-error  of  the  integral  estimate 
V  as 

D  =  (S  —  \L)/N, 
p  =  p  +  D, 

\  =  (N  -2)\/N  +  DDr . 


The  SIR3  can  be  thought  of  as  a  blend  of  a  cubature 
rule  (CR)  and  the  MC  method.  The  SIR3  consists  of 
the  spherical  integration  rule,  which  generates  a  set  of 
2 nx  points  located  on  an  unit  nx  -sphere,  and  the  radial 
integration  rule,  which  governs  spread  of  the  points  within 
the  R"\ 

In  Figure  1,  a  comparison  of  the  points  and  weights 
generated  by  the  CR  and  the  points  and  weights  generated 
by  the  SIR’S  is  shown.  Contrary  to  the  CR,  the  SIR3 
provides  asymptotically  exact  estimates  of  integral  [24], 
i.e.  E[p]  =  p. 

Note  that  the  points  generated  by  the  CR  are  identical 
with  the  points  generated  by  the  UT  with  scaling  parameter 
k  —  0.  In  [13],  it  was  shown  that  the  location  of  points 
generated  by  the  SIR3  when  choosing  a  trivial  orthogonal 
matrix  U  =  I  and  the  parameter  p  as  p2  =  k  +  nx 
matches  the  location  of  points  generated  by  the  UT  which 
uses  the  scaling  parameter  k.  This  fact  lead  to  coining 
the  term  randomized  UT  in  [13]  for  the  SIR3  used  to 
estimate  moments  of  a  random  variable  y  e  R">  related 
to  a  random  variable  x  e  M"r  that  is  given  by  its  mean 
x  and  covariance  matrix  P,  through  a  known  nonlinear 
function  y  =  g(x)  =  [gi(x), . . g„y(x)]r.  The  RUT 
estimates  the  mean  E[y]  =  E[g(x)]  and  the  covariance 
matrix  cov[y]  =  E[(y  —  y)(y  —  $)T  ]  of  y,  and  the  cross¬ 
covariance  matrix  E[(x  —  x)(y  —  y)r], 

IV.  Gaussian  Sum  Global  Filter  with 
Randomized  Unscented  Transform 

The  RUT  introduced  in  the  previous  section  will  be  used 
to  calculate  the  predictive  measurement  moments  (11-13) 


and  the  predictive  state  moments  (21-22)  of  the  Gaussian 
sum  global  filter  described  in  Algorithm  1.  Due  to  paper 
size  limitations,  only  measurement  prediction  moments 
will  be  discussed  in  this  section.  Calculation  of  the  state 
prediction  moments  is  analogous. 

The  simplest  application  of  the  RUT  in  computation 
of  the  predictive  measurement  moments  in  Algorithm  1 
would  comprise  a  SIR3  algorithm  iterative  run  for  each 
of  (V*p_i  means  E^_|[h*(x*)],  A/* |*_i  covariance  matri¬ 
ces  cov^_,[h|t(x*)]  and  Nk\k-]  cross-covariance  matrices 

Pf-f  lL-i-  Naturally,  a  large  number  of  random  parameters 
p  and  random  orthogonal  matrices  Q  would  have  to  be 
drawn  which  would  lead  to  a  substantial  increase  in  com¬ 
putational  costs.  To  increase  the  timeliness  in  the  presence 
of  high  computational  costs,  it  is  possible  to  compute  all 
the  moments  in  parallel.  This  technique  will  be  adopted  in 
the  following  algorithm. 

Additionally,  to  reduce  the  costs  even  further, 
the  algorithm  will  compute  second  raw  moments 
Ejq*_  i  [h*  (x*  )ht  (x*  )T]  and  E^_,[x*h*(x*)T]  instead 
of  the  covariance  and  cross-covariance  matrices 
covKt_1fl,*(x*)]  and  Pj^  jqjfc— l >  respectively.  As  the  last 
step  of  the  algorithm,  the  covariance  and  cross-covariance 
matrices  will  be  obtained  by  a  trivial  combination 
of  the  means  and  second  raw  moments  according  to 
cov(y)  =  E[yyT]  -  E[y][E[y]]T.  This  enables  to  run  a 
parallel  computation  of  all  terms  involved  in  the  Gaussian 
sum  (8). 

For  notational  simplicity,  the  algorithm  will  consider  a 
fixed  number  of  iterations.  Utilization  of  the  error  tolerance 
e  will  be  discussed  later.  Also  for  the  sake  of  clarity,  the 
time  indices  will  be  dropped  in  the  algorithm.  Note  that  p 
and  S  in  Algorithm  2  corresponds  to  E^[  ]  and  E^),JV[  ], 
respectively  in  Algorithm  3.  . 

Algorithm  3:  RUT  for  Gaussian  sum  filter 

Step  1:  Choose  the  total  number  of  iterations  Ntotai  and 
initialize  the  number  of  iterations  as  N  —  0.  Set  the  initial 
values  of  the  measurement  statistics  to  be  calculated: 

.  EU>[h(x)]  =  0n;Xl, 

.  E<2)[h(x)h(x)T]=0„.xnj, 

.  EO)[xh(x)T]  =  0niXnj. 

Calculate  a  square  root  VPU)  of  the  state  predictive 
covariance  matrices  P^ri. 

Step  2:  Repeat  until  N  =  Ntotai,  the  following  loop: 

a)  Set  V  =(V  + 1. 

b)  Generate  a  uniformly  random  orthogonal  matrix 
U  of  dimension  nx  x  nx  and  generate  a  random 
number  p  according  to  p  ~  Chi(nx  +2). 

c)  For  each  of  the  predictive  terms  <V{x;  P(y)} 

do  the  following 

cl)  Compute  a  set  of  points  x;  and  appropriate 
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weights  Wj  according  to 

Xo7)  =  x0),  co0  =  1  -  y_, 

XW=xO')-pU(yPW),-,  COi  =  ^, 
X,(7)+/  -  X0)  +  pU(VPW)/,  (Onx+i  -  COi, 

where  /  =  1,2,...,  n*. 

c2)  Compute  the  values  E*-,^A,[h(x)], 
EO').'v[h(x)h(x)T],  and  E(-/)jA,[xh(x)T] 
of  the  statistics  at  current  iteration 

2  nx 

E0)’,V[h(x)]  =  2>h(X,0}), 

;=0 
2  nx 

EW’>(x)h(x)T]  =  ^fflih(xP)h(xP)T, 

i=0 

2nx 

E0-)^[xh(x)T]  =  X®»-xF)h(xF))T. 

i=0 

Step  3:  Calculate  values  of  the  statistics  as 

N total 

E(;)[h(x)]  =  1  £  Etf).tf[h(x)], 

Af=  1 

N total 

E«[h(x)h(x)T]  =  ^  X  E°)’iV[h(x)h(x)T], 

JV=1 

N total 

E«[xh(x)T]  =  jj  ^  E^’N[xh(x)T]. 

N=  I 

Step  4:  Find  the  covariance  and  cross-covariance  matrices 
as 

COV0)[h(x)]  =  E0)[h(x)h(x)T]  -  E0')|h(x)]E0)[h(x)]T, 
E(7)[(x  -  x(7))(h(x)  -  E0)[h(x)])T]  =  E0)[xh(x)T] 
-x°')E0)[h(x)]T. 


Instead  of  specification  of  a  fixed  total  number  of 
iterations,  it  is  possible  to  monitor  error  of  the  means 
E^’[h(x)],  E(2)[h(x)h(x)T]  and  E^[xh(x)T]  estimates, 
compare  the  error  with  a  pre-specified  threshold,  and  iterate 
until  error  of  all  statistic  estimates  are  lower  than  the 
threshold. 

V.  Numerical  Illustration 


The  measurement  z*  is  influenced  by  the  measurement 
noise  vu  ~  W[o*;  0,  ICC5},  VC  and  the  scalar  parameters 
are  <f> 2  =  0.2  and  ^>3  =  0.5.  The  initial  condition  is  given 
by  a  sum  of  five  Gaussian  pdf’s  p(x 0)  =  S/=i  x 

=AC{(xo;  4J\  PdJ>!  =  Zj=!  0.2  X  ^{xo;  j  -  3,  10}.  The 
predictive  pdf  p(xo|z-1)  is  equal  to  p(x 0). 

For  the  purposes  of  the  Gaussian  sum  global  filters,  we 
calculated  a  three-term  Gaussian  sum  approximation  of 
Ga( 3,2)  distribution  by  means  of  the  expectation  maxi¬ 
mization  algorithm  as 

p(wk)  =  0.29  x  M{wk\  2.14,  0.72} 

+  0.18  x  M{wk\  7.45,  8.05} 

+  0.53  x  M{wk\  4.31, 2.29},  VC 

Performance  of  the  following  state  estimation  methods  was 
compared  in  the  numerical  example: 

•  global  filters: 

-  Gaussian  sum  filter  with  the  RUT  (GSF-RUT), 

-  Gaussian  sum  filter  with  Taylor  expansion  of  the 
nonlinear  equations  [4]  (GSF-TE), 

-  Gaussian  sum  filter  with  unscented  transform  [20] 
(GSF-UT). 

•  local  filters: 

-  EKF, 

-  UKF, 

-  RUKF. 

Within  all  global  filtering  methods,  the  pruning  step  had 
to  be  implemented  to  prevent  exponential  increase  of  the 
number  of  terms  in  the  filtering  Gaussian  sum  pdf.  More 
specifically,  at  each  time  instant,  the  20  highest-weighted 
terms  were  kept  while  the  other  were  discarded. 

The  GSF-RUT  selected  parameters  were  the  maximum 
number  of  iterations  Nmax  —  500  and  the  error  tolerance 
s  —  0.5.  The  efficient  implementation  of  the  RUT  given  by 
Algorithm  3  was  used.  In  this  example,  computational  costs 
of  the  efficient  implementation  amounts  to  approximately 
7%  of  the  computational  costs  of  the  simplest  implemen¬ 
tation  described  in  the  previous  section. 

The  experiments  were  carried  out  using  M  =  1000  MC 
simulations.  Due  to  the  global  property  of  the  estimates 
produced  by  the  GSF-RUT,  five  metrics  were  chosen  for 
comparison  of  the  obtained  results: 

•  Root  Mean-Square  Error  (RMSE)  defined  as 


Consider  the  following  nonlinear  non-Gaussian  system 
[18]  with  one-dimensional  state 


Xk+i  —  4>\Xk  +  1  +  sin(®7rk)  +  Wk  (27) 


with  the  state  noise  Wk  described  by  the  Gamma  pdf 
Ga( 3,2),  Vk,  cj>\  —  0.5,  co  —  0.04  are  scalar  parameters 
and  k  —  1 , . . . ,  60.  The  state  is  observed  by  the  scalar 
measurement  described  by  the  equation 


< h*l  +  Vk, 

faxk  -2  +  0*, 


k  <  30, 
k  >  30. 


(28) 


M 

RMSE*  =  -E  X(%(0  -  x*(/))2, 

N 

where  x*(0  and  x*  ( i )  denote  true  and  estimated  state 
at  the  ;-th  MC  run. 

The  RMSE  metric  provides  an  evaluation  of  the 
estimate  error  expressed  as  the  Euclidean  distance 
between  the  true  state  and  its  estimate.  The  value 
of  the  RMSE  provides  an  absolute  evaluation  of  the 
estimate  error. 
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•  Log  Mean-Square  Error  Ratio  ( log-MSER )  given  by 


M 

log-MSER,  = 

(=1 


tr(Zx-) 

tr(P*|*(/))’ 


where  P*|*(i)  is  the  covariance  matrix  of  the  estimate 
provided  by  the  filter  at  the  i-th  MC  run  and  X,  is 
the  mean-square  error  of  the  estimate.  The  log-MSER 
is  a  scalar  measure  of  the  normalized  non-credibility 
matrix  (NNCM),  NNCM  =  PA“A/2X,P^/2  -  I  [26], 
According  to  the  values  of  the  log-MSER  (posi¬ 
tive/negative)  the  estimator  is  said  to  be  optimistic  or 
pessimistic. 

•  Non-Credibility  Index  (NCI)  defined  as  [11],  [27] 


M 

NCI,  =±  ^l101°gio  ((**('‘)  -Jc,„(i))2/P,|,(i))) 
i=l 

-  101og,o  ((**(!')  -%(i))2/I*)  ]. 

The  NCI  provides  an  evaluation  of  a  relative  esti¬ 
mation  error  [28]  and,  moreover,  it  evaluates  a  self- 
assessment  provided  by  each  filter  in  the  form  of  the 
covariance  matrix  of  the  estimate  error. 

Note  that  NCI  is  a  scalar  measure  dependent  on  the 
distribution  of  the  estimation  error,  while  the  non¬ 
credibility  matrix  is  a  matrix  measure  depending  only 
on  the  actual  and  calculated  second  moments  of  the 
estimation  error  [26]. 

•  Averaged  Normalized  Estimation  Error  Squared 
(ANEES)  [32]  defined  as 


M 

anees,  =  ^  X  ((4  -  4)T(ni*r'(4  -  4))  • 

i=l 

The  ANEES  credibility  measure  has  a  nice  property 
of  dimension  normalization  [27]. 

•  Inaccuracy  defined  as 

K, (p\,P2)=  J  Pi  (X, \Zk )  log 

where  p\  is  the  true  filtering  pdf  and  p2  is  the  filtering 
pdf  obtained  by  the  filters.  The  true  pdf  was  computed 
in  the  form  of  an  empirical  pdf  using  a  particle  filter 
[29]  with  105  samples.  The  global  filters  produced 
the  approximate  filtering  pdf  in  the  Gaussian  sum 
form  (8).  For  the  local  filters,  which  produce  only  the 
filtering  mean  and  covariance  matrix,  the  filtering  pdf 
was  assumed  to  be  Gaussian. 

Considering  the  true  pdf  p\  in  the  empirical  pdf 
form,  the  inaccuracy  [30]  is  a  suitable  measure  of 
discrepancy  between  the  pdf’s  p\  and  p2-  It  equals 
to  the  Kullback-Leibler  distance  between  p\  and  p2 
increased  by  the  Shannon  differential  entropy  of  p\ 
[31]- 

The  RMSE  values  are  given  in  Figure  2,  the  log-MSER 
values  in  Figure  3,  the  NCI  values  in  Figure  4,  the  ANEES 
values  in  Figure  5,  and  the  inaccuracy  in  Figure  6.  Note 


Figure  2.  Time  development  of  the  RMSE. 


Figure  3.  Time  development  of  the  log-MSER. 


that  the  inaccuracy  values  for  the  EKF  and  UKF  were  by 
several  orders  higher  than  the  values  for  the  other  filters 
and,  therefore,  for  clarity  purposes  they  are  omitted  in  the 
figure.  The  inaccuracy  values  for  the  RUKF  ans  GSF-RUT 
were  very  close;  hence  they  are  depicted  on  a  separate 
figure  (Figure  7).  Note  that  the  abrupt  changes  at  k  —  30  are 
caused  by  the  fact  that  at  this  time  instant  the  measurement 
function  becomes  linear  (see  (28)).  As  was  expected, 
the  RMSE  results  show  that  the  GSF-RUT  achieves  the 
smallest  error  of  the  global  filters  and  also  that  global 
estimators  perform  better  than  the  local  ones.  In  terms  of 
the  log-MSER  and  the  NCI,  the  GSF-RUT  provides  the 
best  results  among  the  global  filters.  However,  the  lowest 
values  of  log-MSER  and  NCI  were  achieved  by  the  RUKF, 
which  is  very  surprising  as  it  is  a  local  filter.  The  inaccuracy 
then  again  confirmed  the  expectations  that  global  estimates 
are  closer  to  the  true  filtering  pdf  than  the  local  estimates 
represented  by  a  Gaussian  distribution.  Among  the  global 
estimates  the  results  with  the  lowest  inaccuracy  belong  to 
the  GSF-RUT. 

In  summary,  the  GSF-RUT  achieved  the  best  results  in 
terms  of  the  density  estimates.  However,  this  fact  does  not 
necessarily  imply  that  the  moments  calculated  from  the  best 
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Figure  4.  Time  development  of  the  NCI. 
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Figure  5.  Time  development  of  the  ANEES. 
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Figure  6.  Time  development  of  the  inaccuracy. 


Figure  7.  Time  development  of  the  inaccuracy  for  the  RUKF  and  GSF- 
RUT  (detail). 


EKF 

UKF 

RUKF 

GSF-TE 

GSF-UT 

GSF-RUT 

Time 
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0.4 

0.5 

3.1 

5.2 

7.9 

7.1 

Table  I 

Computational  costs  of  a  time  step  [msec] 


filtering  pdf  estimate  are  also  the  best,  as  was  illustrated 
by  this  example.  Although  the  mean  estimate  achieved 
by  the  GSF-RUT  was,  according  to  the  RMSE,  the  best, 
the  RUKF  accomplished  estimate  error  variance  closer  to 
the  true  one  than  the  GSF-RUT  (measuring  by  the  log- 
MSER).  This  could  also  lead  to  the  claim  that  the  RUKF 
is  a  more  credible  estimator  than  the  GSF-RUT  despite  the 
fact  that  the  GSF-RUT  gives  better  density  estimates.  The 
reason  for  this  situation  may  be  the  fact  that  the  Gaussian 
sum  based  density  estimate  produced  by  the  GSF-RUT  is 
only  an  approximation  of  the  true  one.  For  the  sake  of 
completeness,  the  computational  costs  of  a  time  step  of 
each  algorithm  are  given  in  Table  I. 

VI.  Conclusion 

The  paper  dealt  with  state  estimation  of  nonlinear  non- 
Gaussian  systems  by  Gaussian  sum  filters.  To  achieve 
higher  quality  estimates,  a  Gaussian  sum  filter  was  pro¬ 
posed  where  the  predictive  state  and  measurement  moments 
were  computed  by  the  randomized  unscented  transform 
(RUT).  The  RUT,  which  can  be  seen  as  a  randomized 
version  of  the  unscented  transform,  provides  asymptotically 
exact  estimates  of  the  moments.  To  keep  computational 
costs  of  the  new  filter  low,  an  algorithm  involving  parallel 
calculation  of  the  moments  has  been  proposed  which  is 
termed  the  Gaussian  Sum  Filter  with  Randomized  Un¬ 
scented  Transform  (GSF-RUT).  The  GSF-RUT  filter  was 
illustrated  in  a  numerical  example  and  can  compared  to 
the  local  filters:  extended  Kalman  Filter  (EKF),  Unscented 
KF  (UKF),  and  randomized  UKF  (RUKF)  as  well  the 
global  filters  of  the  Gaussian  Sum  Filter  using  the  Taylor’s 
expansion  (GSF-TE)  and  unscented  transform  (GSF-UT). 
Comparative  result  were  analyzed  in  detail  using  several 


2010 


performance  measures  evaluating  accuracy  of  both  point 
estimates  and  density  estimates.  Three  conclusions  from  the 
comparative  study  were  (1)  in  all  cases,  the  RUT  improved 
the  results  over  the  UT,  (2)  the  RUKF  local  filter  resulted 
in  a  better  log  mean  square  error  and  credibility  index, 
and  (3)  the  global  GSF-RUT  filter  had  the  lowest  absolute 
inaccuracy. 
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